#include <stdlib.h>
#include <stdio.h>
#include <math.h>
/* find the ratio  power spectra
 */
 
int main(int argc, char *argv[]) {
  // USAGE: ratio location filebase #ofseedsofdata
  // Output is written in location
  char filebase[150], fname[200], string[500];
  // open the files in 
  sprintf(filebase, "%s/bias/", argv[1]);
  float col1c[130], col2c[130], avg[130], loc[150], g1[130], ratio[130];
  float col1m[130], col2m[130], col1g[130], col2g[130];
  int i,j, n, ndata;
  FILE *fp;
  
  n = atoi(argv[3]);
  printf("n: %d\n", n);
  for (i=0; i<n; i++) {
    sprintf(fname, "%s%s.mpk.%d", filebase, argv[2], i);
    printf("%s\n", fname);
    fp = fopen(fname, "r");
    ndata=0;
    while (!feof(fp)) {
      if (fgets(string, 499, fp)) {
	sscanf(string, "%g %g %*g", &col1m[ndata], &col2m[ndata]);
	ndata++;
      }
    }
    close(fp);
    
    sprintf(fname, "%s%s.cpk.%d", filebase, argv[2], i);
    fp = fopen(fname, "r");
    ndata=0;
    while (!feof(fp)) {
      if (fgets(string, 499, fp)) {
	sscanf(string, "%g %g %*g", &col1c[ndata], &col2c[ndata]);
	ndata++;
      }
    }
    close (fp);
    
    sprintf(fname, "%s%s.gpk.%d", filebase, argv[2], i);
    fp = fopen(fname, "r");
    ndata=0;
    while (!feof(fp)) {
      if (fgets(string, 499, fp)) {
	sscanf(string, "%g %g %*g", &col1g[ndata], &col2g[ndata]);
	ndata++;
      }
    }
    close (fp);
    // now take the ratio et al and write the output
    sprintf(fname, "%s%s.ratiocm.%d", filebase, argv[2], i);
    fp = fopen(fname, "w");
    for (j=0; j<ndata; j++) {
      if (col2m[j]>0) {
	ratio[j]=col2c[j]/col2m[j];
        avg[j]=(col1c[j]+col1m[j])/2;
        fprintf(fp, "%g %g\n", avg[j], ratio[j]);
      }
    }
    close (fp);
    
    sprintf(fname, "%s%s.ratiogm.%d", filebase, argv[2], i);
    fp = fopen(fname, "w");
    for (j=0; j<ndata; j++) {
      if (col2m[j]>0) {
	ratio[j]=col2g[j]/col2m[j];
	avg[j]=(col1g[j]+col1m[j])/2;
	fprintf(fp, "%g %g\n", avg[j], ratio[j]);
      }
    }
    close (fp);
    
    sprintf(fname, "%s%s.stoch.%d", filebase, argv[2],i);
    fp = fopen(fname, "w");
    for (j=0; j<ndata; j++) {
      if (col2m[j]>0) {
	ratio[j]=col2c[j]*col2c[j]/(col2g[j]*col2m[j]);
	avg[j]=(col1g[j]+col1m[j]+col1c[j])/3;
	fprintf(fp, "%g %g\n", avg[j], ratio[j]);
      }
    }
    close (fp);
  }
}
